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Abstract 

When the Vlasov equation is investigated numerically using the method of 
test particles, the particle-particle interactions that inevitably arise in the 
simulation (but are not present in the Vlasov equation itself) result in an 
accumulation of errors which eventually drive the collection of test particles 
toward a state of classical thermal equilibrium. We estimate the rate at which 
these errors accumulate. 


The Vlasov equation plays a central role in classical (and semiclassical) time-dependent 
mean field theory, and has been used to model a wide variety of many-body processes, 
from the gravitational iV-body problem |I|, to plasma physics @], to nuclear dynamics 0. 
While the content of the Vlasov equation is conceptually simple — interactions among 
many particles are replaced by a common mean-field potential — solutions are harder to 
come by, and must in general be sought numerically. This is often accomplished with the test 
particle method: a swarm of numerical particles is used to simulate a distribution /(r,p,t) 
in one-body phase space, and the mean-held potential in which these test particles evolve is 
obtained from this distribution. Thus, while the Vlasov equation replaces a physical many- 
body problem with the self-consistent evolution of a one-body phase-space distribution, 
the test particle method in turn replaces the Vlasov equation with a numerical many-body 
problem. This raises the issue of convergence: for a given number of test particles, and 
over a given length of time, how closely can we expect the evolution of /(r, p, t) as obtained 
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by the test particle method, to resemble the true evolution under the Vlasov equation? In 
the limit of arbitarily long evolution time, we can certainly expect the test particle method 
to fail. In that limit, the unavoidable interactions between individual test particles will 
drive the swarm of numerical particles toward a Boltzmann distribution of energies, whereas 
under the Vlasov equation there are no particle-particle interactions, and / typically does 
not evolve toward classical thermal equilibrium. 

A relevant example of this disagreement arises in the application of the Vlasov equation to 
nuclear dynamics, where the Pauli principle is imposed by insisting that, initially, / < 4/h^ 
everywhere in phase space. Under the Vlasov equation, this condition is preserved exactly 
with time; with the test particle method, however, classical thermalization occurs, and the 
Pauli condition is violated. (See Fig. 4 of Ref. [Q, I, for an illustration.) 

The eventual thermalization of test particles may occur on a time scale longer than 
that in which one is interested. Nevertheless, it is indicative of a general process, whereby 
interactions between the test particles introduce errors which drive the numerical solution 
away from the actual solution of the Vlasov equation. It it thus important to obtain an 
estimate of the rate at which these errors accumulate. Such an estimate is the goal of the 
present brief note. 

We will restrict ourselves mainly to self-consistent potentials which are local functions of 
particle density, although a brief discussion of long-range forces will be presented at the end. 
Working with a simple schematic model, we will argue that errors in the particle energies 
accumulate diffusively, and we will solve for the functional dependence of the associated 
diffusion constant De, in terms of physical and numerical parameters. For the specihc case 
where Gaussian smoothing is used to obtain the particle density p from the positions of the 
test particles, we obtain a more quantitative prediction for De- Finally, we compare our 
theoretical predictions with numerical results. 

The Vlasov equation is explicitly given by 

^ + {/,i/} = 0, (1) 
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(2) 


where denotes the ordinary Poisson bracket, and the Hamiltonian H is 

The notation Hj(r) is meant to indicate that the mean-held potential U{r) is a functional of 
the one-body phase space distribution /. Often (e.g. when the physical interactions between 
particles are independent of momentum) the functional dependence of H on / reduces to a 
dependence only on the density p in ordinary space: 


(7/(r) ^ (7,(r) 

(3) 

P(r) = J dp f(r,p)- 

(4) 


Throughout this paper, we will assume, for simplicity, that this is the case. Note that if 
U is linear in p, then it may be expressed in terms of a two-body physical interaction V 12 . 
f/p(r) = / (ir'p(r') Vi 2 (r, r'). In general, however, U need not be linear in p, i.e. the mean 
held need not arise from two-body interactions. 

The implementation of the test particle method involves two tasks: (1) evolving each of 
the N test particles in the presence of the time-dependent potential U ; and (2) constructing 
U from the positions of the particles at any instant in time. The hrst is straightforward, 
involving simply the numerical integration of Hamilton’s equations of motion. The second 
task requires the particle density p(r,t), which is obtained by smearing the position of each 
point particle with a localized folding function g: 

p(r,t) = — '^g(r-ri(t)). (5) 

i=i 

Here A is the number of physical particles, whereas the sum runs over the test particles, 
located at positions rj(t) at time t. g{r) is a function localized in a volume around the 
origin, and normalized to unity: / dr g{r) = 1. The parameter a thus measures the distance 
over which we smear out the particle positions. Gaussian folding functions are commonly 
used. 

To estimate the rate of accumulation of errors introduced by interactions among the 
test particles, let us consider a simple model in which our many-particle system is conhned 
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within a box of volume V. Furthermore, let us take the functional dependence of U on p to 
be local: 

[/p(r) = [/(p(r)). (6) 

That is, the potential at r depends only on the density of particles at that point; this 
corresponds to zero-range interactions among particles, and is commonly used to model 
short-range interactions such as nuclear forces (see e.g. the Skyrme parametrization |^). It 
is important to distinguish here between the mathematical problem one is trying to solve 
(propagation under the Vlasov equation), and the numerical method used to solve it: even 
if the potential f/p(r) which enters into the Vlasov equation is exactly local — as indeed we 
are assuming — the interactions between test particles in a numerical implementation will 
necessarily have finite range, due to the smearing which is employed to extract a smooth 
density p(r) from the positions of a finite number of test particles. 

Now, consider an initial phase space distribution /o(r, p) corresponding to an ensemble 
of monoenergetic particles distributed uniformly throughout the box, with an isotropic dis¬ 
tribution of momenta. Explicitly, this has the form /o(r, p) oc S(p—Po)Ob(j'), where p = |p|, 
and ©^(r) is equal to 1 (0) if r is inside (outside) the box. As can be seen by inspection, 
this phase space distribution is a stationary solution of the Vlasov equation, thus under the 
Vlasov equation the ensemble of particles remains exactly monoenergetic. Our strategy now 
will be to investigate how such an initial distribution evolves under a numerical simulation 
using test particles. Specifically, after a time At, what is the amount AB by which the 
energy of a typical test particle has strayed from its initial value? The growth of AB with 
At will then be a measure of the accumulation of error inherent in the test particle method. 

Let us take our N test particles — all given the same initial speed v = pQ/m — to be 
distributed randomly throughout the container, and choose a so that V/N <^V. This 
will result in a reasonably smooth numerical density p(r, t), without smearing over too large 
a volume of the box. We can express this density as 

p(r,t) = po + Sp{v,t), (7) 
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where po = is the physical density which we are trying to simulate, and (5p(r, t) rep¬ 
resents the fluctuations around po due to the flnite number of test particles. To gauge the 
typical size of 6 p, note that the value of p at a given point is roughly equal to {A/N) n/a^, 
where n is the number of test particles within a volume of the point in question. (The 
factor A/N is a conversion factor between the density of test particles and the density of 
physical particles.) On average, n will be given by no = Na^/V, with fluctuations of size 
around this average. These considerations yield the following expression for the typical 
size of the fluctuations Sp: 



It should be clear as well that Sp{ri,t) and 6 p{r 2 ,t) will be correlated only if ri and r 2 
are within a distance ~ a of one another. Furthermore, at a given location r, the value of 
Sp{r, t) will be correlated over a time tc ~ cr/v, since that is a typical time over which a test 
particle remains within a volume element of r. 

Thus, our numerical density p(r, t) fluctuates in space and time around an average value 
Po = AjV^ where the size of the fluctuations is given by Sprms ~ Po /and these fluc¬ 
tuations are correlated over a distance a, and a time tc ~ a/v. Let us now make use of 
this picture to determine what happens to a given test particle evolving under the potential 
U{p{r,t)) computed from this numerical density. 

We first expand the potential U around its value at po: 


U{p{jA)) = f/(po+ 5p(r,t)) ^ Uo + UgSpirA), (9) 


where Uq = U{po) and U/ = ^(po)- Thus, like the numerical density p, the potential U 
fluctuates in space and time around an average value {Uq). The typical rate at which U is 
changing, at a fixed point r, is determined by the typical rate of change of Sp: 

kc/(r,«)Uc/i^. (10) 

Now consider a single test particle i moving under this time-dependent potential. From 
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Hamilton’s equations, the rate of change of the total energy Ei of the particle is exactly the 
value of dU/dt along its trajectory: 

E,(t) = (11) 

This function Ei{t) is autocorrelated over a time scale tc ^ cr/v, which is considerably shorter 
than a characteristic time scale associated with the particle’s motion in the box (e.g. the 
traversal time across the length of the box). The change in energy AEi is thus the time 
integral of a function Ei{t) which fluctuates rapidly, with short time correlations; this implies 
that AEi evolves diffusively. The associated diffusion constant De is then the time integral 
of the auto-correlation function of Ei{t). Approximating this integral by the product of the 
mean-square value of Ei{t) with the correlation time tc, we have, using Eq.pT], 


Df 


E, 


2 cr 




2 V 

a 


Finally, using 6prms ~ Po/^/Eq, and uq = Na^/V, we get 


( 12 ) 


De ~ . h . 1. (13) 

iV a 

Thus after a time At S> we can expect the energy of our test particle to have changed 
by an amount 


AE ~ {DEAt)E'^, 


(14) 


with De given by Eq.|^ above. Eqs.^ and |T^ together describe the accumulation of error 
in the energy of a typical test particle, and thus constitute our main result. 

Note that we have written De sls the product of three factors, the first of which contains 
only physical quantities, while the other two depend on purely numerical parameters: the 
smearing parameter a, and the number of test particles per physical particle, N/A. The 
prediction that the error accumulates more slowly for larger values of N/A is expected; this is 
the benefit of using more test particies. Eq.|^ predicts that one gains even more by increasing 
the vaiue of the smearing parameter a: De oc As pointed out by Reinhard and Suraud 
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0 (II), this should not come as a surprise: a larger smearing effectively suppresses the 
interaction between individual test particles, thus slowing the rate at which energy gets 
exchanged. Of course, smearing distorts the mean held itself; therefore too much of it, while 
suppressing errors due to test particle interactions, will result in an inaccurate simulation of 
the Vlasov equation. Ultimately, one wants a large enough so that AE remains small over 
the time scale of physical interest, but not so large as to distort the inhomogeneities that 
are physically present in the mean held. 

As mentioned earlier, one expects that in the long run the test particles thermalize. This 
ought to happen on a time scale r over which each test particle has had the opportunity to 
change its energy by an amount comparable to the average particle energy, mv‘^/2. Thus, 

~ (15) 


Combining this with Eq.p!3|, we obtain for the thermalization time scale 

4 N 

T ~ - • (J • — 

(t/avo ^ 


(16) 


In numerical experiments aimed at studying the relaxation toward thermal equilibrium under 
the test particle method, Reinhard and Suraud have found that doubling the value of a 
“gains more than an order of magnitude in the relaxation time” (Ref. |^, p. 227); this is in 
agreement with our prediction here that r scales like Furthermore, these authors have 
predicted that r cx N/A, and have conhrmed this numerically. 

A few comments are now in order. First, in a realistic test particle simulation, the 
particles are held together by the mean held itself, rather than being artihcially conhned 
within a box. Nevertheless, the mechanism by which the test particles exchange energy with 
one another remains the same, therefore the result derived within the context of our simple 
model ought to hold in the more realistic situation as well. 

Next, while our main result predicts how the growth of errors scales with the various 
parameters involved, a more quantitative estimate will depend on the details of how the 
test particles interact with one another. For instance, the use of a gaussian folding function. 
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g{r) = (27rcr^)“^/^ exp(—r^/2(j^), allows for an explicit evaluation of Sprms- This leads to an 
expression for De which has the form of Eq.[T^, but with a numerical factor in front. 

Alternatively, for gaussian folding functions one can evaluate (within the linear approxi¬ 
mation) the amount of energy exchanged in a given collision between two test particles, in 
terms of impact parameter. The further assumption that different particle-particle collisions 
are uncorrelated leads (after some work) to a diffusion coefficient 




(17) 


It is encouraging that this approach, which differs somewhat from that leading to Eq.p!B|, 
nevertheless yields the same functional dependence of De on the various quantities involved. 

Finally, our assumption that the dependence of 17 on p is local (Eq.|) was made both for 
the sake of simplicity, and because our original motivation to study this problem arose from 
the application of the Vlasov equation to nuclear dynamics, where short-range physical forces 
lead to a local Up. However, in many physical applications of the Vlasov equation one deals 
with long-range forces (e.g. Coulombic and gravitational forces), therefore it may be useful 
to extend the analysis of the present work, to include non-local mean-held potentials. It is 
interesting in this context to note that Chandrasekhar [Q has made a detailed calculation of 
the time scale Te required for binary stellar interactions to drive a self-gravitating system 
of many stars (e.g. a galaxy) toward thermal equilibrium. His result, translated into our 
notation, is Te ~ /{Gm?y‘p, where we have removed dimensionless factors.^ Now, the 

gravitational potential at a given point in the galaxy is roughly U —NGm^/R, where R 
is a distance scale characterizing the size of the galaxy, and N is the number of stars. If, for 
purposes of comparison with Eq.|T^ above, we write U' UIp U/NR then we get 


2 2 
m V 

{U'fp 


■ R'. 


(18) 


^ including the logarithm of a quantity which is essentially the ratio of kinetic to potential energy 
for a typical star in the system. 
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This has the form of our Eq.|T^, only with a replaced by R, and without the factor N/A 
(since there are no “test particles”). The fact that the size of the entire galaxy, R, appears in 
place of our smoothing parameter a, suggests that — in comparison with short-range forces 
— long-range forces such as gravity strongly suppress the collisional relaxation toward ther¬ 
mal equilibrium. This is consistent with numerical Endings: in semiclassical simulations of 
nuclear dynamics, the (undesirable) approach to classical Boltzmann equilibrium often takes 
place on a time scale comparable with the mean-field dynamics in which one is interested 
IP. In contrast, simulations of the many-body gravitational problem evolve rapidly to a 
near-static “collisionless equilibrium”, which differs from a microcanonical distribution [0 . 

We now present the results of numerical experiments which we have performed to test our 
predictions. We simulated the schematic model discussed above — a gas of particles confined 
to a box — where the box was taken to be a cube of volume 1000 (in arbitrary units) with 
periodic boundary conditions, and the mean field potential was taken as U{p) = —Rp + Ap^. 
A = 200 physical particles of mass m = 1 were assumed, the initial speed of each particle 
was set to n = 1, and a gaussian folding function was used. In each simulation, we allowed 
the gas of particles to evolve for a time At = 2.5, and we followed the growth of the mean 
square change in test particle energy, ((AE)^) = {1/N)Y,f=i{AEiy‘^ over this time. For 
each simulation, we found ((AF)^) to grow linearly with timeQ, in agreement with Eq.[l^. 
To extract a numerical energy diffusion coefficient De, we divided the final value of ((AE)^) 
by At. 

We ran two sets of simulations. In the first set, the smearing parameter was held fixed at 
CT = 1, and the number of test particles was varied from N = 2000 to iV = 10000. In Fig.|l| 
we plot the resulting values of as a function of N/A on a log-log plot. From Eq.|T^, we 
expect the data to fall along a straight line of slope -1; the solid line gives the best fit of the 
data to a line of this slope. In the second set of simulations, the number of test particles 


^aside from the inevitable short quadratic growth at the start 
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was held fixed at 4000, and the smearing parameter was varied from cr = .5 to cr = 1. Fig.|^ 
plots the resulting values of as a function of a on a log-log plot, and the solid line gives 
the best fit of the data to a line of slope -4. In both figures we find good agreement between 
the prediction of Eq.|^ and the numerical results. If we allow the slope of the lines to vary 
as well, then a best fit of straight lines to the two sets of data yields slopes of —1.08 and 
—4.54, instead of —1 and —4. Thus, the discrepancy between the predicted slope and the 
best fit to numerical data is on the order of 10% in each case. 

Since we used a gaussian folding function in our numerical simulations, it is interesting 
to compare the results directly to the quantitative prediction of Eq.|^. That prediction is 
depicted by the dotted line in each of the figures. We see that Eq.0 overestimates the value 
of He by a factor of nearly 2. We believe that this discrepancy is due to our neglect of 
correlations between different collisions. (Since motion in our mean field is highly regular, 
we can expect that the energy exchanged in a given collision between two test particles 
will be somewhat correlated with the energy exchanged at the next collision of those two 
particles.) 

To conclude, we have argued in this paper that when the Vlasov equation is studied 
numerically using the method of test particles, the inevitable interactions between test par¬ 
ticles introduce errors which accumulated diffusively. Our main result — a prediction of 
how the associated diffusion constant scales with the available numerical parameters — was 
shown to agree well with the results of computer simulations. 

We acknowledge conversations with P.G. Reinhard and J. Randrup which stimulated and 
improved this work. This work was supported by the Department of Energy under Grant 
No. DE-FG06-90ER40561. 
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FIGURES 

FIG. 1. The energy diffusion constant De as a function of N/A, the number of test particles 
per physical particle. The heavy dots show the values of De extracted from numerical experiments, 
the solid line shows the best fit of these points to a straight line of slope -1, and the dotted line 
shows the quantitative prediction of Eq.^. 

FIG. 2. De as a function of smearing parameter a. The solid line is a best fit of the data to a 
straight line of slope -4; the dotted line represents EqJT7|. 
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